On the energy spectrum of strong magnetohydrodynamic turbulence 
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The energy spectrum of magnetohydrodynamic turbulence attracts interest due to its fundamental 
importance and its relevance for interpreting astrophysical data. Here we present measurements of 
the energy spectra from a series of high-resolution direct numerical simulations of MHD turbulence 
with a strong guide field and for increasing Reynolds number. The presented simulations, with 
numerical resolutions up to 2048^ mesh points and statistics accumulated over 30 to 150 eddy 
turnover times, constitute, to the best of our knowledge, the largest statistical sample of steady 
state MHD turbulence to date. We study both the balanced case, where the energies associated 
with Alfven modes propagating in opposite directions along the guide field, E'^{k±) and E~{k±), are 
equal, and the imbalanced case where the energies are different. In the balanced case, we find that the 
energy spectrum converges to a power law with exponent —3/2 as the Reynolds number is increased, 
consistent with phenomenological models that include scale-dependent dynamic alignment. For the 
imbalanced case, with E^ > E~ , the simulations show that E~ oc k^^^^ for all Reynolds numbers 
considered, while E'^ has a slightly steeper spectrum at small Re. As the Reynolds number increases, 
E^ flattens. Since E^ are pinned at the dissipation scale and anchored at the driving scales, we 
postulate that at sufficiently high Re the spectra will become parallel in the inertial range and scale 



as -E oc -E oc fc' 



-3/2 



Questions regarding the universality of the spectrum and the value of the 



"Kolmogorov constant" are discussed. 



I. INTRODUCTION 

Astrophysical plasmas are typically magnetized and 
turbulent, with turbulent fluctuations spanning a 
tremendous range of scales in which the energy spec- 
trum follows a power law scaling [e.g., [l|, Incom- 
pressible magnetohydrodynamics (MHD) provides the 
simplest theoretical framework for studying magnetized 
plasma turbulence. The precise form of the MHD tur- 
bulence spectrum is crucial for a variety of processes in 
astrophysical systems with extended inertial intervals, 
such as plasma heating and wave-particle interactions, 
which are sensitive to small variations in the spatial scal- 
ing of the fluctuations [e.g., l^-lsj. The incompressible 
MHD equations take the form 



T • V z 



V-z± = 



(1 



where = v ± b are the Elsasser variables, v is the 
fluctuating plasma velocity, b is the fluctuating magnetic 
field (in units of the Alfven velocity), Va = 'Bo/^/A^Tpo is 
the Alfven velocity based upon the uniform background 
magnetic field Bq, P = {p/po + ^^/2), P is the plasma 
pressure, po is the background plasma density, i/ is the 
fiuid viscosity (which, for simplicity, we have taken to 
be equal to the magnetic diffusivity) and represent 
forces that drive the turbulence at large scales. It can be 
shown that in the limit of small amplitude fiuctuations, 
and in the absence of forcing and dissipation, the system 
describes non-interacting linear Alfven waves with dis- 
persion relation w^(k) = ik^^VA- The incompressibility 



condition requires that these waves be transverse. Typi- 
cally they are decomposed into shear Alfven waves (with 
polarizations perpendicular to both Bq and to the wave- 
vector k) and pseudo- Alfven waves (with polarizations in 
the plane of Bq and k and perpendicular to k). 

Nonlinear interactions (or collisions) between counter- 
propagating Alfven wave packets distort the packets, 
splitting them into smaller ones until a scale is reached 
when their energy is converted into heat by dissipa- 
tion 0. The efficiency of the nonlinear interaction is 
controlled by the relative size of the linear and nonlin- 
ear terms in equation ([T]): The regime in which the lin- 
ear terms dominate over the nonlinear terms is known as 
weak MHD turbulence, otherwise the turbulence is called 
strong. The Fourier energy spectrum of MHD turbulence 
can be derived analytically only in the limit of weak tur- 
bulence [e.g., However, it has been demonstrated 
both analytically and numerically that the energy cas- 
cade occurs predominantly in the plane perpendicular 
to the guiding magnetic field [1, [13, [H, [Tl|, which en- 
sures that even if the turbulence is weak at large scales 
it encounters the strong regime as the cascade proceeds 
to smaller scales. Although weak turbulence may exist 
in some astrophysical systems [e.g..[l3l. [T7l - [T9| . magnetic 
turbulence in nature is typically strong, for which an ex- 
act analytic treatment is not available. In this case, high- 
resolution, well-optimized numerical simulations play a 
significant role in guiding our understanding of the tur- 
bulent dynamics. This provides the motivation for the 
present work. 

The ideal MHD system conserves the Elsasser energies 
E+ ^ \ J{z+)^d^x and = \ J{z-y(fx (equivalcntly, 
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the total energy E = E'" + = \ j {v'^ + b'^)(Px = E+ + 
E~ and the cross-hcUcity — J[v ■ h)(fix — £'+ — E— 
are conserved). The energies E'^ and E~ cascade m a 
turbulent state toward small scales due to the nonlinear 
interactions of oppositely moving z"*" and z~ Alfven pack- 
ets. MHD turbulence is called balanced when the energies 
carried by oppositely moving fluctuations E^ are equal, 
and it is called imbalanced when they arc not the same. 
MHD turbulence in nature and in the laboratory is typi- 
cally imbalanced. For instance, this is the case when the 
turbulence is generated by spatially localized sources, as 
is the case in the solar wind where more Alfven waves 
propagate away from the Sun than towards it. The inde- 
pendent conservation of the two Elsasser energies (com- 
pared to only one conserved energy in hydrodynamics) 
has a profound consequence for the MHD dynamics [e.g., 

Mm- 

In this work we present the results of a series of direct 
numerical simulations of MHD and Reduced MHD for 
balanced and moderately imbalanced turbulence and in- 
vestigate how the scalings of the Elsasser spectra behave 
as the Reynolds number is increased. We also present the 
first high-resolution direct comparison of simulations of 
MHD vs RMHD turbulence, demonstrating that the lat- 
ter model completely captures the turbulence dynamics 
of strong MHD turbulence at roughly half the computa- 
tional cost of a full MHD simulation. 

This paper is organized as follows. In section |ll] we 
briefly describe the most recent phenomcnological efforts 
to understand scaling laws in MHD turbulence, particu- 
larly in the imbalanced case. In section IIIII we describe 
the numerical set up and the parameter regime for our 
simulations. In section II VI we show measurements of the 
energy spectrum from a series of numerical simulations 
with varying Reynolds numbers. In section |V] we show 
measurements of scale-dependent dynamic alignment and 
establish its relation with the —3/2 scaling of the energy 
spectrum. In section |Vl] we discuss the approach to the 
universal regime and the universality of Kolmogorov's 
constant in MHD. We show that dynamic alignment in- 
troduces a new robust scale-dependent quantity that en- 
ters the definition of the energy spectrum and uniquely 
sets the Kolmogorov constant. We propose that this new 
quantity is a consequence of cross-helicity conservation. 
Finally, in section IVIII we discuss our results. 

II. MHD TURBULENCE PHENOMENOLOGY 

For strong MHD turbulence, Goldreich & Sridhar [s^l 
argued that the pseudo- Alfven modes are dynamically 
irrelevant for the turbulent cascade (since strong MHD 
turbulence is dominated by fluctuations with k± ^ 
the polarization of the pseudo- Alfven fluctuations is al- 
most parallel to the guide field and they are therefore 
coupled only to field-parallel gradients, which are small 
since fcy <^ k±). If one filters out the pseudo- Alfven 
modes by setting zj}^ = 0, it can be shown that the re- 



sulting system is equivalent to the Reduced MHD model: 

(^^ T • V||^ z± + (zT • Vi) z± = -ViF 

+iyV^z^+if, (2) 

We note that in RMHD the fluctuating fields have only 
two vector components, but that each depends on all 
three spatial coordinates. Moreover, because the are 
assumed incompressible ( V • = 0) , each field has only 
one degree of freedom which is more commonly expressed 
in terms of stream functions in the more standard form 
of the RMHD equations [Ullll. 

Conservation of both the Elsasser energies means that 
once an imbalance has been created it cannot be de- 
stroyed by the MHD dynamics. It is also well known that 
decaying MHD turbulence, affected only by the dissipa- 
tion, becomes increasingly more imbalanced with time 
[e.g., [13, [H, [13 ■ Several analytic and numerical studies 
have shown that imbalance is also an inherent property of 
driven MHD turbulence even if the turbulence is forced 
without introducing a net imbalance at the largest scales 
- the turbulent domain spontaneously fragments into lo- 
cal imbalanced domains where the cross helicity is either 
positive or negative [2l| - |26l . [28l . [29| . 

In imbalanced domains, the directions of the magnetic 
and velocity fluctuations are not independent, rather, 
they are either aligned or counter-aligned to a certain 
degree [s^. The organization of such a domain is the 
following: the directions of both the magnetic and ve- 
locity fluctuations vary within a small angle (comparable 
to the alignment angle) throughout the domain, while 
their amplitudes change predominantly in the direction 
normal to their polarizations. Such positively and nega- 
tively aligned domains appear to be the building blocks 
of MHD turbulence, whether it is balanced overall or not. 

The origin of such domains can be qualitatively under- 
stood from the conservation of energy and cross-helicity 
in an ideal MHD system. When small dissipation is 
present and the system is unforced, it can be argued that 
energy decays faster than cross-helicity. This selective 
decay would eventually lead to Alfvcnization of the flow, 
that is, to progressively stronger alignment (or counter- 
alignment, depending on the initial state) between the 
directions of the magnetic and velocity fluctuations, e.g., 
[20I [21I m, nil . In a perfectly aligned (counter-aligned) 
state either z"*" or z~ is identically zero, and the nonlinear 
interaction vanishes. In a driven state, characterized by 
strong nonlinear interaction and a constant energy flux 
over scales, the alignment cannot be perfect. Rather, 
it turns out that alignment depends on the scale, the 
smaller the scales the better the alignment. Below we will 
demonstrate this phenomenon in numerical simulations. 
From a more qualitative point of view, one can argue 
that whenever a partly aligned domain appears, nonlin- 
ear interaction inside such a domain gets reduced, and 
its evolution time increases compared to non-aligned do- 
mains. Therefore aligned domains persist longer, which 
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explains the tendeney of a turbulent flow to exhibit such 
self-organization. These aligned domains are the domains 
where the essential energy of the turbulence is contained, 
and they are typically well seen in numerical simulations. 
Solar wind observations also show that globally balanced 
turbulence is made up of locally imbalanced patches at 
all scales 

In the aligned or imbalanced domains, the Elsasser en- 
ergies are unequal, and one can ask whether their spec- 
tra have to be the same. This raises questions of whether 
MHD turbulence is universal and scale-invariant. Indeed, 
if imbalanced domains have different spectra that depend 
on the degree of imbalance, their superposition may not 
have a universal scaling. 

Phcnomenological treatment of strong imbalanced 
MHD turbulence is complicated by the fact that one can 
formally construct two time scales for the nonlinear en- 
ergy transfer: The times of nonlinear deformation of the 
packets at some spatial scale A arc ^ -^/^Ji which 
can be significantly different in the case of strong imbal- 
ance, [e.g., [10, [131 ■ In recent years, several phcnomeno- 
logical models attempting to accommodate this differ- 
ence have been proposed. However, the theories have 
generated conflicting predictions because they use differ- 
ent assumptions regarding the physics of the nonlinear 
energy cascade. For example, the theory by Lithwick 
et al. (ssl concludes that in the imbalanced regions the 
Elsasser spectra have the scalings E'^{k±) cx E~{k±) cx 

kj^^^; the same spectra were also suggested by Beres- 
nyak and Lazarian [s^. The theory by Chandran 
proposes that the spectra of E^{k±) and E^{k±) are 
different depending of the degree of imbalance, while the 
theories by Perez and Boldyrev [1^ and Podesta and 
Bhattacharjee [4l[ find that the spectra of E^{k±) and 
E~{k±) have different amplitudes but the same scalings 

E+{kj_) (X E-{k±) (X k-^/^. 

One would expect that numerical simulations could 
clarify the picture. However, the first numerical sim- 
ulations of strongly imbalanced MHD turbulence [e.g., 
[H, [H, \^ also produced conflicting results regarding 
which power law E^ should follow. The conflicting nu- 
merical findings apparently reflect the fact that imbal- 
anced MHD simulations require significantly more com- 
putational effort compared to the balanced cases (43| . 
This happens since in the imbalanced domains the non- 
linear interaction is depleted and the Reynolds and mag- 
netic Reynolds numbers are reduced. This can be for- 
mally seen from the fact that, in a strongly imbalanced 
domain with z"*" ^ z~ , the z'^ field is advected by a low- 
amplitude z~ field, and therefore z~^ becomes directly 
affected by the dissipation at smaller wave- vectors (com- 
pared with the balanced case), which reduces its incrtial 
interval. Now, z~ is advected by a strong but z^ 
is significantly affected by the dissipation, so the inertial 
interval of z~ becomes spoiled as well. 

In order to produce large inertial intervals simultane- 
ously for both Elsasser fields when strongly imbalanced 
domains are present in the flow, one therefore needs to 



have a significantly higher Reynolds number as com- 
pared to the balanced case. However, as one increases 
the Reynolds number, one needs to increase the numeri- 
cal resolution in order to appropriately resolve the small 
scales and to make sure the numerical run is stable. 
Therefore, the larger the imbalance, the larger the numer- 
ical resolution required to describe correctly the Elsasser 
spectra. Fortunately, it has been argued that Reduced 
MHD can be used to investigate the universal properties 
of MHD turbulence, which offers the advantage that an 
RMHD simulation can be achieved at half the cost of an 
MHD simulation. 



III. NUMERICAL SETUP 

We solve the MHD equations ^ and their RMHD 
counterpart ((2) in a periodic, rectangular domain with 
aspect ratio x Ly, where the subscripts denote the 
directions perpendicular and parallel to Bq, respectively. 
We set = 27r, L\\/Lj_ = 6 or 10 and Bq = Se^. A 
fully dealiased 3D pseudo-spectral algorithm is used to 
perform the spatial discretization on a grid with a res- 
olution of N'^ X iV|[ mesh points. We note that the do- 
main is elongated in the direction of the guide field in 
order to accommodate the elongated wave-packets and 
to enable us to drive the turbulence in the strong regime 
while maintaining an inertial range that is as extended 
as possible (see [44| ) . This is a physical requirement that 
should be satisfied no matter what model system, full 
MHD or reduced MHD, is used for simulations. 

In the case of reduced MHD though, when the zjj^ 
components are explicitly removed, the resulting system 
(21) is invariant with respect to simultaneous rescaling 
of the background field Bq and the field-parallel spatial 
dimension of the system, if one neglects the dissipation 
terms. Therefore, for any strength of the background 
field Bq ^ I, one can rescale the field to = 1 and the 
field-parallel box size to Ly = L±, that is, conduct the 
simulations in a cubic box. We should note however that 
the dissipation terms in ^ are not invariant and they 
should be changed accordingly under such rescaling. 

To save on computational cost we have reduced the 
field-parallel numerical resolution for some simulations, 
i.e., the numerical grid is anisotropic with L\\/N\\ > 
L±/N_i_. This is appropriate since the energy cascade 
proceeds much faster in the field-perpendicular direc- 
tion, and the energy spectra decline relatively slowly in 
the field-perpendicular direction and relatively fast in the 
field-parallel direction. Energies at large fcy are therefore 
reduced and a lower field-parallel resolution is not ex- 
pected to alter the behavior of the spectra in the incrtial 
interval. An isotropic resolution with the value imposed 
by the field-perpendicular dynamics would therefore be 
wasteful. 

We should however caution that a reduced resolution 
(or, equivalently, unreasonably high Reynolds number 
for a given resolution) may contaminate the dissipative 
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physics, even if the inertial interval is unaffected. For 
example, if the precise scaling behavior in the dissipa- 
tion interval is of interest, as is the case for extended 
scaling laws such as the dynamic alignment angle, some- 
what smaller Reynolds numbers may need to be chosen. 
As a general rule, whether the numerical simulations are 
conducted to investigate the inertial or the dissipation in- 
terval, a resolution study must be performed in order to 
establish the optimal Reynolds number for a given task. 
In particular, it has to be verified that increasing the nu- 
merical resolution while keeping the physical parameters 
such as Reynolds number, forcing mechanism, etc. un- 
changed does not affect the studied spectra, e.g., (45| . 
This point will be illustrated below in the balanced case. 

The turbulence is driven at the largest scales by collid- 
ing Alfven modes [s^. We drive both Elsasser popula- 
tions by applying statistically independent random forces 
f"'' and f~ in Fourier space at wave- numbers 2'k/L± < 
k± < 2{2tt/L±), /cji = 27r/L||. The forces have no com- 
ponent along z and are solcnoidal in the xy-planc. All of 
the Fourier coefficients outside the above range of wave- 
numbers are zero and inside that range are Gaussian ran- 
dom numbers with amplitudes chosen so that Vrms ^ 1- 
The individual random values are refreshed indepen- 
dently on average approximately 10 times per turnover of 
the large-scale eddies. The variances aj. = (|f±P) control 
the average rates of energy injection into the z"*" and z~ 
fields. We take a'^ > a~ and in the statistically steady 
state we measure the degree of imbalance through the pa- 
rameter h= {E+-E-)/{E+ + E-) ^ /E. Thus h = 
corresponds to balanced turbulence and h ^ \ defines 
maximally imbalanced turbulence. Time is normalized 
to the large scale eddy turnover time tq = L±/ (2-KVrms)- 
The field-perpendicular Reynolds number is defined as 
Re± = Vrms{Li_l2'iT)/v 1/v. In order to accommodate 
the reduced field-parallel resolution we have also modi- 
fied the diffusion operator in equations ((Ij and ©, i.e., 
we have replaced with v{dxx + dyy) + v\\dzz- 

The system is evolved until a stationary state is 
reached, which is confirmed by observing the time evo- 
lution of the total energy of the fluctuations. The data 
are then sampled in intervals of the order of the eddy 
turnover time. All results presented correspond to av- 
erages over 30-150 samples for each run. As shown in 
Table [ml we conduct a number of MHD and RMHD 
simulations in the balanced and imbalanced regime in 
order to investigate the scaling of the energy spectra as 
the field-perpendicular Reynolds number increases. 



TABLE I: Simulation Parameters. 



Case 


Regime 




iV|| 


h 


L\\/L± 


Re±_ 




RBla 


RMHD 


512 


256 





6 


2400 


V 


RBlb 


RMHD 


512 


512 





6 


2400 


V 


RBlc 


RMHD 


512 


512 





6 


1800 


V 


RB2a 


RMHD 


1024 


256 





6 


6000 


2.5i/ 


RB2b 


RMHD 


1024 


1024 





6 


6000 


V 


RB2c 


RMHD 


1024 


1024 





6 


3200 


V 


RB2d 


RMHD 


1024 


1024 





6 


1800 


V 


RB3a 


RMHD 


2048 


512 





6 


15000 


2.5i/ 


RB3b 


RMHD 


2048 


2048 





6 


15000 


V 


RB3c 


RMHD 


2048 


2048 





6 


9000 


V 


RB3d 


RMHD 


2048 


2048 





6 


5700 


V 


RIl 


RMHD 


512 


256 


0.45 


10 


2200 


V 


RI2 


RMHD 


1024 


256 


0.5 


10 


5600 




RI3 


RMHD 


2048 


512 


0.5 


10 


14000 


2.5u 


MBl 


MHD 


512 


256 





10 


2200 


V 


MB2 


MHD 


1024 


256 





10 


5600 




Mil 


MHD 


512 


256 


0.5 


10 


2200 


V 


MI2 


MHD 


1024 


256 


0.5 


10 


5600 


2.5i/ 


MI3 


MHD 


2048 


512 


0.5 


10 


14000 


2.5!^ 



TABLE II: Summary of the numerical runs with different nu- 
merical resolutions and different Reynolds numbers. There 
is no particular scheme used to choose the Reynolds number 
for a given resolution other than to ensure that the studied 
scaling properties are well demonstrated and the numerical 
runs are stable. 



1.0 
















0.1 




' 2048^x512 

1024^x256 

512^x256 


\ \ ■ 









1 10 100 



FIG. 1: Total field-perpendicular energy spectrum in bal- 
anced RMHD as the Reynolds number increases (Cases RBla, 
RB2a, RBSa). 



IV. MEASUREMENTS OF THE ENERGY 
SPECTRUM 

The field-perpendicular energy spectrum is obtained 
by averaging the angle-integrated Fourier spectrum, 

E{k^) = 0.5(|v(ki)|2)fci + 0.5(|b(ki)|2)fc^, (3) 



over field-perpendicular planes in all samples. Identify- 
ing the inertial range in numerical simulations with lim- 
ited resolution is generally difficult, due to the relatively 
modest separation between the forcing and dissipation 
scales that current super-computers can afford. For in- 
stance, a measurement of the turbulence spectrum for 
a single Reynolds number is not enough to ensure that 
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FIG. 2: Total field-perpendicular energy spectrum in bal- 
anced MHD as the Reynolds number increases (Cases MBl, 
MB2). 
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FIG. 3: Energy spectra E'^{k±) and E {k±) in imbalanced 
RMHD as the Reynolds number increases (Cases RIl, RI2, 
RI3). 
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FIG. 4: Energy spectra E'^{k±) and E {k±) in imbalanced 
MHD as the Reynolds number increases (Cases Mil, MI2, 
MIS). 



[T] to |4] the inertial range is identified by the flat regions 
of the spectra compensated by k^^^, which extend fur- 
ther to the right with increasing Reynolds number (and 
resolution) . 

Figures [T] and [5] show the total field-perpendicular en- 
ergy spectrum E{k±) in the balanced regime for the 
RMHD and MHD cases, respectively. The RMHD and 
MHD spectra are remarkably similar, confirming that 
the pseudo-Alfven modes are dynamically insignificant 
and that the RMHD approximation is valid. In both 
cases the total energy spectrum remains of the form 
E(k±) kj^^'^ as the Reynolds number increases, with 
the inertial range starting at /c « 4 and extending up to 
A: > 30 in the highest Reynolds number case. In neither 
RMHD or MHD is there any evidence of a build up of 
energy close to the dissipative wave-numbers-often re- 
ferred to as a bottleneck effect- with both spectra falling 
off smoothly in the dissipative range. 

Figures [3] and [1] show the field-perpendicular Elsasser 
spectra in the imbalanced regime for the RMHD and 
MHD cases, respectively. Again the behavior of both 
spectra in the RMHD and MHD regimes are very sim- 
ilar. In both cases, it is seen that while E~ keeps the 
scaling E~(k±) ~ k^^'^ as the Reynolds number in- 
creases, the scaling of E'^{k±) is more difficult to pin 
down. Indeed, both the RMHD and MHD results for 
Re = 2200 yield a steeper spectrum for E'^{k±)^ with 
an exponent possibly nearer to —5/3 than —3/2. How- 
ever, we believe that there is no real significance to the 
value of —5/3 here, the exponent is simply steeper than 
—3/2. Indeed, in both cases, as the Reynolds number is 
increased E~^(k±) appears to flatten, which means that 
£'^(fc_L) lias not fully established the universal scaling be- 
havior yet. Since E^{kj_) and E^{k±) are pinned (i.e., 
converge to each other) at the dissipation scales and are 
anchored (i.e., independent of the Reynolds number) at 
the driving scales, we postulate that at sufficiently high 
Re (where the inertial range is extensive) the spectra will 
become parallel in the inertial range and attain the scal- 



ing E^{kj_) - k' 



-3/2 



Numerical tests of this prediction 



must await a significant increase in computational power. 



V. MEASUREMENTS OF DYNAMIC 
ALIGNMENT 

An important test that can be performed in the pre- 
sented simulations concerns the so-called dynamic align- 
ment angle. This angle is defined by the following ratio 
of the two specially constructed structure functions [2^ : 



the simulated turbulence has converged to the asymp- 
totic universal scaling. Instead, one carries out a set 
of numerical simulations with increasing resolution and 
Reynolds number. The spectra are then compensated by 
the different phcnomenological predictions and the pre- 
ferred model is distinguished by the best fit. In Figures 



9{l) = 



(|^v^(I)x^b^(l)|) 
(|<5v^(l)||Jbx(l)|) ^ 



(4) 



where <5vj_(I) and (5bx(l) are the field-perpendicular ve- 
locity and magnetic field increments, respectively, cor- 
responding to the field-perpendicular scale separation 1. 
(We note that in definition ^ we have assumed that 
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the angle is small, and hence no distinction between 6{l) 
and sm9{l) is made. Hereafter, by 6(1) we will always 
understand the quantity (jl}). 

As proposed in [23l. [23| the alignment angle 9{l) has a 
nontrivial scaling with I, which may explain the observed 
—3/2 scaling exponent of the energy spectrum. As dis- 
covered in [45j . the scale dependent dynamic alignment 
exists not only in the inertial interval, but it also extends 
into the dissipation range and is limited only by the grid 
size of the numerical scheme. We will demonstrate that 
the alignment angle scaling provides a sensitive test prob- 
ing the turbulent cascade deep in the dissipation interval. 
In particular we will see that if the simulated dissipation 
range is under- resolved (e.g., as a result of the use of 
too large a Reynolds number or strongly anisotropic res- 
olution), the dynamic alignment can be easily spoiled at 
the dissipation scales even if it is present in the inertial 
interval. 

The measurements of the alignment angle are pre- 
sented in Figure [Sj The first panel shows three simu- 
lations (RB2d,c,b in Table IIII[) performed at the same 
numerical resolution of 1024'^ but with different Reynolds 
numbers Re 1800, 3200, 6000. Plots for Re = 
1800, 3200 show a remarkable property of the alignment 
scaling: It extends deep down into the dissipation region, 
practically up to the scale of the numerical discretization, 
independently of the Reynolds number (see also [ist). 
However, this behavior is spoiled if the Reynolds number 
is pushed to very high values, at which the dissipation 
interval becomes under-resolved. In this case, the scaling 
starts to degrade at large wave-numbers, as is seen in the 
case Re = 6000. 

The alignment scaling is however restored back to its 
original value if the numerical resolution is increased to 
2048^^, so that the dissipation scales become well resolved 
again. This is seen from comparison of the plot for RB2b 
(1024^, Re = 6000) in the first panel of Figure [5] with 
the plot for RB3d (2048^, Re = 5700) in the second 
panel. Further increase of the Reynolds number in the 
second panel of this figure demonstrates that the align- 
ment scaling is stable up to Re = 9000 (RB3c, 2048^), 
however, it starts to degrade at large wave-numbers for 
higher Reynolds numbers Re = 15000 (RB3b, 2048^), in 
complete analogy with the behavior depicted in the up- 
per panel of Fig. [5] at smaller resolution. The alignment 
angle is spoiled even more in the run RB3a (2048^ x 512, 
Re = 15000) in the same figure where we simultaneously 
decrease the field-parallel numerical resolution, making 
the dissipation interval even more under-resolved. Note 
however, that in both the first and the second panels of 
Fig. [5l the heaviest distortion of the alignment behavior 
occurs in the dissipation region, while the inertial inter- 
val (approximately contained between the two vertical 
lines) is relatively unaffected. This may explain why an 
under-resolved dissipation interval is not manifest in the 
scaling of energy spectra, as seen in Figure |6l 

Fig. [7] shows three well resolved simulations with nu- 
merical resolutions increasing from 512'^ to 1024"^ to 




0.2 1 . . 

0.001 0.010 0.100 

l/2n 




0.2 L...^ ... I 

0.001 0.010 0.100 
l/2n 

FIG. 5: Measurements of the dynamic alignment angle (|4]) 
vs scale I in balanced RMHD. Upper panel: simulations 
RB2d (solid), RB2c (dashed), RB2b (dash-dotted) on 1024^ 
mesh points. Lower panel: simulations RB3a (dash-triple- 
dotted) on 2048^ X 512 mesh points, RB3b (dash-dotted), 
RB3c (dashed), RB3d (solid) on 2048^ mesh points at differ- 
ent Reynolds numbers. The dynamic alignment scaling ex- 
tends well into the dissipation range, up to scales close to 
the grid cell (roughly I ~ 3 grid cells). When the Reynolds 
number is pushed to very high values (so that the dissipation 
interval becomes under-resolved) or the numerical resolution 
in the field-parallel direction is reduced, the alignment-angle 
scaling degrades at small scales. The vertical lines show the 
approximate boundaries of the inertial interval (cf. Fig. ((8|). 
The straight dotted line has a slope of 1/4. 



2048^^. We observe that the scaling interval of the align- 
ment angle becomes progressively longer and its scaling 
index stays close to the predicted value 1/4 with 
little or practically no dependence on the Reynolds num- 
ber This means that we observe a truly universal 
scaling behavior of the dynamic alignment. The lower 
panel of Fig. [7] shows the same curves where the spatial 
scale is normalized by the dissipation length. We observe 
that the flattened parts of the curves at small scales do 
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FIG. 6: Energy spectra for runs RB3a (solid) and RB3d 
(dash). Simulation RB3a on 2048^ x 512 has an unresolved 
dissipation at the expense of longer inertial interval. Simula- 
tion RB3d is performed at lower Re to capture alignment in 
the dissipation region, with a shorter inertial range. 



not overlap under such rescaling, which supports our ob- 
servation mentioned above that the extent of the scaling 
interval is not defined solely by the dissipation scale, but 
rather depends on the numerical discretization step. 



VI. ENERGY SPECTRUM: KOLMOGOROV 
CONSTANT AND DISSIPATION SCALE 

For a more complete study of the energy spectrum, 
one can also evaluate the amplitude of the spectrum and 
the dissipation scale for each simulation and verify that 
they agree with a given phenomenology. Since our spec- 
tral scaling conforms to the phenomenology of Boldyrev 
[23I [23 | , we now study in more detail the scaling associ- 
ated with this model. First, we need to derive the ex- 
pression for the energy spectrum, which is done in the 
following way [2^. The time of nonlinear interaction at 
field-perpendicular scale A in this model \s t ^ X/ {v\6\), 
where vx denotes the typical (rms) velocity fluctuations, 

= Oo{X/L±y^^ is the scale-dependent alignment angle 
between magnetic and velocity fluctuations, which was 
studied in the previous section, and 6q is the typical align- 
ment angle at the outer scale (forcing scale) L±. The rate 
of energy cascade is then evaluated as e = v^9\/X, from 
which it follows that Eik^) - e2/3(0^/2.y^)-2/3/c-^/^ 
One however notices that the amplitude of the energy 
spectrum is not uniquely defined in this equation, since 
the outer-scale quantities and Lj_ essentially depend 
on the forcing routine. This is understood from the fol- 
lowing example. Assume that the large-scale force drives 
only unidirectional Alfvcn waves z"*", for which v is per- 
fectly aligned with b and = 0. Then the wave energy 
will grow without bound, since the nonlinear interaction 
leading to the energy cascade and eventual dissipation at 
small scales is absent. 

Even when a particular forcing routine is specified, the 
definitions of the values of 9q and L± are still subjective 
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FIG. 7: Measurements of the dynamic alignment angle (|4]) 
in balanced RMHD. The frames show numerical simulations 
with increasing Reynolds number and numerical resolution 
with properly resolved dissipation ranges (runs RBlc (dash- 
dotted), RB2c (dashed), and RB3d (solid)). In the lower 
plot, the scale I is rescaled by the dissipation length (see sec- 
tion |Vl] for precise definitions). It can be seen from here that 
the region of scale dependent dynamic alignment increases as 
smaller scales are made available by increased numerical res- 
olution. The extent of the alignment region is not limited 
by the dissipation scale, but rather depends on the grid size 
of the numerical scheme. The straight dotted line has the 
slope 1/4. 



since they essentially rely on the outer-scale properties 
of turbulence rather than on the measurements of the 
inertial interval. We now propose that this problem can 
be remedied in an efficient way. For that we notice that 
there exists a well-defined quantity that is remarkably 
stable (scale-independent) in the inertial interval: 



(5) 



where 6{l) is defined in see the discussion in the 
preceding section. In this definition one can use any 
scale I from the inertial interval or dissipation interval if 
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the numerical simulations are well resolved. A somewhat 
simpler rule can be used in numerical (or observational) 
studies, where one docs not have to know a priori what 
scales correspond to the inertial interval and does not 
have the luxury of having the plot in Fig. [7] available. In 
this case / in formula (O can be chosen to be the Taylor 
micro-scale based on either the magnetic or the velocity 

fluctuations, / = frms/|V X v|rms or / = &rms/|Vxb|rms, 

assuming the magnetic Prandtl number is of order one. 
We therefore propose the following normalization of the 
energy spectrum: 

where A is defined by ([5]). The scale A that is defined 
solely through the inertial-interval quantities, incorpo- 
rates the essential information about the cross-helical 
structure of MHD turbulence. It is not uniquely defined 
by the outher scale of the turbulence, rather it also de- 
pends on the large-scale driving mechanism. Therefore, 
the inertial-interval energy spectrum is defined by the two 
quantities e and A. characterizing the energy cascade rate 
and the level of cross-helical organization of the fiow. The 
presence of the two quantities characterizing the spec- 
trum of MHD turbulence (as oppose to only one quantity 
in liydrodynamic turbulence) is the manifestation of the 
two conserved quantities cascading toward small scales 
in MHD turbulence: energy and cross-helicity. 

We expect that the constant Ck in ^ may be "univer- 
sal," that is, largely independent of the character of the 
driving, analogous to the Kolmogorov constant in hydro- 
dynamical turbulence. This constant can be measured in 
our simulations in the following way. First, we specify 
I that we use to measure the alignment scale A in ([5|). 
According to our plots in Figs. [5] and [71 we may choose 
I = 0.07L±, say, as a scale belonging to the inertial in- 
terval and not yet affected by the numerical resolution 
effects. Then, for simulations RBla, RB2a, RB3a we 
find A = 1.34L_L, 1A1L±, 1A8L±, respectively. 

The dissipation rate can be evaluated based on the 
energy spectrum ^ as follows: 
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E{kii,kj_){iykl^ + i'\\kf)dk\\dk^ 



(7) 



Our numerical results confirm that the integral of 
leads to a negligible correction to the dissipation rate, 
and therefore it can be omitted, and we can use the 
field-perpendicular spectrum E(k±) = J E{k\\,k±)dk\\. 
Then, for simulations RBla, RB2a, RB3a we find: e = 
0.15,0.15,0.16. 

The dissipation scale can be found (or defined) based 
on the energy spectrum. Omitting the dimensionless con- 
stants, we then accept, by definition, 

^-2/9^1/9^2/3^ 



V 



(8) 



We can demonstrate that our simulations agree with 
this scaling by plotting the energy spectra in the bal- 
anced case (RBla, 2a, 3a) versus the wave- vector normal- 
ized with the dissipation scale (|S]), where we measure 
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FIG. 8: Upper panel: Total field-perpendicular energy spec- 
trum E{k±) in balanced RMHD turbulence for different 
Reynolds numbers (Cases RBla, RB2a, RB3a). The wave- 
number is normalized by the dissipation scale ([8} and the en- 
ergy is compensated by A~^''^e~^''^k^/^ . The rescaled curves 
collapse onto each other (up to the forcing scale) revealing 
the universal functional form of the energy spectrum. Lower 
panel: The scaling of the length of the inertial interval with 
the Reynolds number. Good agreement with the phenomeno- 
logical model ((6] [8]) is observed. 



the dissipation rate directly from the simulations via (O , 
and the alignment scale from ([5]). The top frame in fig- 
ure [8] shows that in this case the dissipative region starts 
around krj « 0.1, independent of the Reynolds number. 
The extent of the inertial range, defined as the ratio be- 
tween the scale Iq at the beginning of the inertial range 
(from figure [U k ~ 4 and hence Iq ~ L±/8) and the dis- 
sipation scale Id ~ Lj_/{2kd) — 5r]L±^ where k^ = 0-1 /rj 
from figure H]) , increases up to one decade in the RB3a 
case 



61|. 



Note that with the wave vector normalized 
with the single parameter rj^ the whole spectra collapse 
onto each other, thus providing additional evidence that 
the universal functional behavior of the spectrum is ob- 
tained in our simulations. The lower plot in figure [8] 
shows that the length of the inertial range increases as 
lo/^d ^ Re^^^, also in good agreement with the estimate 
for the dissipation scale (|8]). The "Kolmogorov constant" 
Ck can be evaluated from the upper plot as Ck 2. 



VII. DISCUSSION 

We have presented results from state-of-the-art di- 
rect numerical simulations of balanced and imbalanced 
driven MHD turbulence. The simulations are achieved 
at the extremely large numerical resolution of 2048^ 
and the longest running time, with many runs span- 



9 



ning more than a hundred eddy turnover times in the 
steady state. The simulations were performed using two 
pseudo-spectral codes, one solving the MHD equations 
and the other solving the RMHD equations. In the- 
ories and simulations of MHD turbulence, it has long 
been argued that RMHD provides a correct and accu- 
rate framework for investigating the universal properties 
of MHD turbulence both in the weak and strong turbu- 
lence regimes. We have presented a direct comparison of 
high-resolution numerical simulations of MHD vs RMHD 
turbulence using two independently developed pseudo- 
spectral codes with identical parameters. It is shown 
that in the strong turbulence regime, in both the bal- 
anced and imbalanced state, the energy spectrum of the 
Elsasser variables in MHD and RMHD are in remarkable 
agreement (for details of a lower resolution comparison, 
including the individual velocity and magnetic spectra 
and the alignment angle, see |4^). These results are of 
essential value for MHD turbulence research, as simulat- 
ing MHD turbulence can be accomplished using RMHD 
codes that generally incur a smaller computational cost. 
In the balanced case, the simulated energy spectra of 
and E~ show a clearly identifiable inertial range, 
consistent with a slope of kj^^^ for both and E~ . It 
is observed from Figures [1] and [5] that the compensated 
energy spectra show a flat region that extends as the 
Reynolds number is increased. This is consistent with 
previous, lower resolution simulations of strongly mag- 
netized MHD turbulence, e.g., fl6ll47l45^. In the imbal- 
anced case, the interpretation of the numerical results is 
not as straightforward. Figures [3] and 0] show that the en- 

3/2 

ergy spectrum of E~ remains reasonably close to fcj^ , 
only slightly changing its overall amplitude for small 
Reynolds numbers. As for the spectrum, the com- 
pensated spectrum shows a slope slightly steeper than 
—3/2 which however flattens as the Reynolds number in- 
creases. Another observation from the large Reynolds 
number imbalanced numerical simulations is that the 
spectra of E'^ and E~ are "anchored" at large scales 
and "pinned" at the dissipation scale. From these re- 
sults we propose that the energy spectra of E~^ becomes 
asymptotically closer to k'^^'^ as the Reynolds number 
is increased. Much higher resolutions, exceeding the ca- 
pabilities of today's supercomputers, are required to con- 
clusively demonstrate this conjecture. 

Finally, during the refereeing process, our attention 
was drawn to rece nt p ublications by the group of Beres- 
nyak & Lazarian [53l - l5a |. in which the authors address 
issues similar to the ones contained in this paper. Most 
of the conclusions of those papers appear to be at odds 
with our s ( and with similar results or other groups, 
G-g-, O l47H5l| ). We however note that the actual nu- 
merical results presented in [ssl - ls^ agree with ours in 
the range of scales that we study, while they differ from 
ours at very large wavenumbers, e.g., fc > 50 in the runs 
with highest resolution. Beresnyak & Lazarian suggest 
that the true inertial interval exists only at these large 
wavenumbers where they perform their measurements 



of the scaling relations. The formal cause of the dis- 
agreement of our conclusions with those by Beresnyak & 
Lazarian is thus the numerical measurements being per- 
formed in essentially different regions of the phase space. 

The question however remains as to what causes the re- 
sults of the numerical simulations by Beresnyak & Lazar- 
ian to disagree with ours at small, sub-inertial scales. Ac- 
cording to our analysis, the answer is the following: the 
fc-space intervals on which references f53l - [55| base their 
conclusions are significantly affected by numerical effects 
due to the numerical setup they use. It is not appropriate 
in their simulations to use those intervals for addressing 
cither the inertial or the dissipation regimes. We how- 
ever note that the dissipation-range dynamics and the 
behavior of the numerical solution of the MHD equations 
close to the numerical cutoff is an interesting and not well 
studied question. It is therefore worth addressing the dif- 
ferences between our simulations and those by [53l455l | in 
more detail. Since such analysis is not the main objective 
of the present work, we have presented the corresponding 
discussion in the Appendix. 
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APPENDIX: Numerical study of MHD tubrulence 
at subrange scales 

In this Appendix we comment on the numerical re- 
construction of solution of the MHD equations at small 
scales, that is, scales within the dissipation range and 
close to the numerical cutoff in fc-space (the dealiasing 
cutoff in a pseudo-spectral code) . Recent publications 
by Beresnyak & Lazarian [53l - l55| found that the en- 
ergy spectrum in this region (roughly corresponding to 
/c > 50 in their highest-resolution runs) can have a pecu- 
liar structure that is inconsistent with the structure and 
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the scahng found in our numerical simulations. Much 
confusion was created by the suggestion by [sj, [HHl that 
this high-A: region is, in fact, the true inertial interval of 
MHD turbulence, while the region that is studied in our 
works (corresponding to 4 < /c < 30 in highest-resolution 
runs) is a "non-converged" forcing-dominated region. 

It is therefore useful to address the small-scale numer- 
ical solution of the MHD equations in more detail and 
to relate our findings and conclusions to those presented 
in the recent works by Beresnyak and Lazarian j53l - l55l |. 
These references claim that the energy spectral index of 
MHD turbulence is —5/3 and that there is no conclusive 
evidence for dynamic alignment in the numerical results. 
In discussing what could lead to such (in our opinion, 
erroneous) conclusions it is useful to distinguish two fac- 
tors. One is related to differences that arise because the 
simulations by Beresnyak & Lazarian that allegedly are 
identical to ours, in fact are not identical at all because 
of differences in the details of the numerical setup. The 
other is related to the methods that are used to analyze 
the results and, ultimately, support one claim or another. 
Both play a role in the origin of the disagreement. 

First, we concentrate on issues that result from the 
different setup of the numerical simulations. In our pre- 
vious publications (e.g., p3. 46]) we have discussed at 
length those aspects of the simulation design that are es- 
sential for accurately capturing the physics of the strong 
turbulent cascade. It is not necessary to repeat those 
discussions here, however, it is important to point out 
that many of the simulations of Beresnyak and Lazarian 
[53l - l55j differ from ours through their choice of numeri- 
cal hyper-dissipation, significantly smaller viscosities for 
a given numerical resolution, and a considerably smaller 
statistical ensemble from which averages are computed. 
Each of these factors is potentially detrimental for the 
observation of the correct scaling behavior. For example, 
the measurements of the alignment angle that are shown 
in Figure 3 of reference |54| and Figure 2 of reference 
[ssj lead Beresnyak to conclude that dynamic alignment 
is not present in MHD turbulence as the alignment angle 
saturates, i.e. flattens as a function of / at small /, when 
the Reynolds number increases. However, those plots 
exhibit a behavior that is qualitatively similar to that 
displayed in our Figure EJ where insufficient numerical 
resolution is demonstrated to affect the alignment angle 
at small scales. It therefore reasonable to conclude that 
the observed flattenin g of the alignment angle in the sim- 
ulations of references |54l . Issj is an artifact of unresolved 
dissipation scales and, possibly, part of the inertial-range 
scales, rather than a physical effect. 

The influence of hyper-dissipation may be similarly as- 
sessed from comparing the energy spectra obtained in 
our work with the energy spectra obtained in, say, ref- 
erence [13] . Our spectra in Figures [TJ [51 [5] & H] exhibit 
an extended interval with the scaling identifled as 

the inertial interval, followed by a steep decline, identi- 
fied as the dissipation range. The spectra in Figure 2 
of reference [541 also show an extended interval with the 



scaling k~^/'^ (interpreted in reference [13] as a "non- 
converged" range) followed at large wavenumbers by a 
very short steepening (interpreted in [H^] as the "inertial 
interval") and then flattening and ultimate cut-off. In 
our view, such an interpretation is incorrect; the spec- 
tral behavior observed in reference [13] close to the dis- 
sipation region is not a property of the inertial inter- 
val, but rather is evidence of the so-called bottleneck ef- 
fect that is expected when numerical hyper-dissipation 
is present. Indeed, as discussed in references [Ha, Ist} . 
an energy spectrum abruptly terminated in fc-space by 
hyper-dissipation or by other Galerkin-type truncation 
mechanisms, exhibits an inertial interval followed by a 
pseudo-dissipation region (steepening of the spectrum), 
then by a partly thermalized region (a rise in the spec- 
trum), and then by a far dissipation range (ultimate cut- 
off). The measurements presented in References [s^. [55| 
are consistent with such spectral behavior, which moti- 
vates a natural explanation of their results as an inertial 
interval with the —3/2 scaling, modified by a substantial 
bottleneck effect close to the dissipation scales. More- 
over, a thermalization brought about by sharp termina- 
tion of the spectrum in the fc-space tends to decorrelate 
small-scale fiuctuations, which otherwise would remain 
strongly aligned throughout the dissipation interval, cf. 
our Figure [7| This is also consistent with the significant 
loss of dynamic alignment at small scales that is observed 
in references [13, [lal ■ 

A more detailed comparison of our results can be made 
with those MHD simulations by Beresnyak [Hj that em- 
ploy a physical Laplacian dissipation (simulations R8 & 
R9 in [55[). By evaluating the Reynolds numbers for 
those calculations in the same way that it is done in our 
work. Re ~ VrmsL/{2'!iv) with Vrms ~ 1, we find that 
simulation R8, with a resolution 768^ mesh points, is per- 
formed at the Reynolds number Re ~ 8000, while calcu- 
lation R9 (resolution 1536^) is performed at Re « 20000. 
According to our results in Figure [5l in the simulations 
with a resolution of 1024^ mesh points the dissipation 
interval is under- resolved already at Re « 6000, while in 
the 2048'^ simulations the dissipation interval is under- 
resolved at Re « 15000. Thus, the runs R8 k R9 of refer- 
ence [ssj that are most similar to ours have lower numer- 
ical resolutions while higher Reynolds numbers. There- 
fore, they have essentially unresolved dissipation intervals 
and, possibly, parts of the inertial intervals. 

The lack of resolution at the bottom of the inertial 
intervals in the simulations R8 &: R9 can also be seen 
from the alignment- angle curves shown in Figure 2 of 
[ssj . Under the rescaling applied in that figure, the curves 
should approach each other in the inertial interval, as 
they do in our Fig. [71 lower panel. In contrast, one can 
see only a short region of in Figure 2 of [111 (runs R8 & 
R9) where the curves approach each other, approximately 
within the range 20 < l/v 40. Apparently, this is the 
only piece of the inertial interval that is resolved, and 
in this interval the scaling exponent of the angle indeed 
approaches 1/4, see Figure 3 in reference [Hj, as expected 



11 



according to our results. 

We now turn to the second factor that contributes 
to the differing conclusions drawn by the Beresnyak & 
Lazarian group, namely the method of analysis. We re- 
call that the objective is to determine the scaling be- 
haviour within the inertial range. Concerning the energy 
spectrum, we assess whether the numerical data preferen- 
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tially supports E{k±) oc kj_^^^ or E{k±) cx kj°''^ directly 
by compensating the numerical data by k^^'^ and by /c^/"^ 
in turn. For the correct model, the inertial range then 
corresponds to the range of scales over which the com- 

3/2 

pcnsatcd spectrum is flat. We always find that k^ 
provides the better fit, with the inertial range starting 
at fc_L 4 and extending up to k± > 30 at highest res- 
olution. As the Reynolds number increases, numerical 
convergence is demonstrated by the fact that this region 
maintains its amplitude and scaling and increases in ex- 
tent to larger wave numbers, see, e.g. our Figure [TJ 

In contrast, Beresnyak [55| uses an indirect method 
to select the preferred spectral exponent. He uses the 
two phenomenological models that describe the inertial 
range characteristics to predict the dissipation scales (ry), 
plots the compensated spectrum as a function of the di- 
mensionless wavcnumbcr krj, and identifies the preferred 
model as that which displays the better convergence 
properties at large wavenumbers krj as the Reynolds num- 
ber increases. Figure 1 of Beresnyak [55| lead him to con- 
clude that it is the —5/3 model that displays the better 
convergence properties at large k. 

It can be shown, however, that the convergence at 
small scales observed in [s^, is a simple artifact of 
the numerical setup adopted in [sj, [H^, rather than a 
physical effect. To explain this, we note that any discrete 
numerical scheme solves only the corresponding discrete 
algebraic equations. If the numerical setup is done cor- 
rectly, the numerical solution approximates the physical 
one independently of the discretization step. If, however, 
a special numerical setup is adopted where rj is rigidly 
tied to the grid size such that r/N is kept the same in 
all runs (as is done in [sj, [s^ ) , then the numerical so- 
lution plotted as a function of kr] is always affected by 
the discretization in the same way, thus consistently re- 
producing the same small-scale numerical effects that are 
present in the setup. The convergence at large kr] is then 
the convergence among solutions of the given numerical 
scheme, which should not be confused with the conver- 
gence to the physical solution. 

To illustrate this effect in our simulations we replot 
the spectra presented in Fig. |8] choosing the Kolmogorov 



-5/3 



normalization scale r]K4i 



.3/4,-1/4 



Due to a partic- 



ular choice of viscosities in our runs depicted in Fig. [8l 
in this case rjKU happens to double every time the res- 
olution decreases by a factor of 2, thus ensuring that 
Vk4i^± = const, see Fig. [3] It is therefore not surprising 
that all the curves converge in the vicinity of the numeri- 
cal dealiasing cutoff corresponding to krjx^i ~ 0.8, while 
they do not converge in the inertial interval and in the 
most of the dissipation interval. A similar, by design. 




FIG. 9: The spectra presented in Fig. [8] rescaled with a new 
parameter 7^x41 = i^^'^'^e"^''*. For this particular choise, 77^41 
is proportional to the step of numerical discretization, that is, 
N±'qK4,i ~ const. According to our estimate of the onset of 
the dissipation region in Fig. [S] the corresponding region in 
the present plot starts at k±riK4,i ~ 0.04, while the numerical 
dialiasing cutoff that is always imposed at k± = N±/3, is seen 
at k±riK4,i ~ 0.8. We observe that the convergence is present 
in the visinity of the dealiasing cutoff, k±riK4,i ^ 0.3, while it 
is absent in the inertial interval and in most of the dissipation 
interval. 
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FIG. 10; Comparison of the spectra obtained in numeri- 
cal runs RB2c and RB3a, rescaled by the Kolmogorov dis- 
sipation scale 77x41 = !^^^*e~^^* (upper panel) and by the 
dissipation scale obtained in the dynamic alignment theory 
77DA = !/^''^A^/^e~^/^, see (jSj (lower panel). In these runs, 
the numerical discretization is not related to the dissipation 
scale as N±ri = const. With no spurious numerical conver- 
gence imposed by the numerical setup, the —5/3 model does 
not fit the data, while the dynamic alignment model shows a 
good agreement in the inertial interval and in the dissipation 
range, up to the scales where the numerical efi'ects eventually 
become significant. 



convergence is present in Fig. 2 of [54| and Fig. 1 of [5f 
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Such convergence at very small scales is a spurious nu- 
merical effect, which does not reflect the convergence of 
the physical solutions, and which cannot give preference 
to any phenomenological model. When the viscosities 
in different runs do not conform to the special condition 
riNj_ = const, the spurious convergence disappears, and 
the —5/3 model does not fit the data, while the —3/2 
model still provides a good fit in the incrtial and dissipa- 
tion intervals, see Fig. [TO] 

We therefore conclude that the numerical simulations 
by Beresnyak & Lazarian group [ssl - lssj are likely signifi- 



cantly affected by numerical effects at small scales where 
their measurements are performed. This is notwithstand- 
ing the statements made in [ssl - fssi that the simulations 
are resolved in those works. These statements, in our 
opinion, are not supported by the factual numerical data 
presented in these papers. Until the effects of hyper- 
dissipation are better understood and numerical conver- 
gence is demonstrated in the settings of [ssi - lssj , it is hard 
to assess fully the degree to which the numerical findings 
of [53l - l55| can be compared with our results, or with sim- 
ilar results of other groups [e.g.. HgI. I47l - l49l [5]|. 
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